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Abstract 

Among Monte Carlo techniques, the importance sampling requires fine 
tuning of a proposal distribution, which is now fluently resolved through it- 
erative schemes. The Adaptive Multiple Importance Sampling (AMIS) of 
Cornuet et al. (2012) provides a significant improvement in stability and Ef- 
fective Sample Size due to the introduction of a recycling procedure. How- 
ever, the consistency of the AMIS estimator remains largely open. In this 
work we prove the convergence of the AMIS, at a cost of a slight modifi- 
cation in the learning process. First numerical experiments exhibit that this 
modification might even improve the original scheme. 
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1 Introduction 

The aim of Monte Carlo techniques is to approximate a target distribution Il(-) on 
some space with a weighted sample. Namely, they output a system of particles 
x e e X (indexed by e e E), with their weights co e e [0;oo). Then, the discrete 
measure 

n £ (-) = YjOj e d x x-) 
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serves as an approximation of the target II(-). And the Monte Carlo scheme is 
said to be consistent if, for a large class of functions if/ : — » R, the sum 
J i^(x)n £ (djc) = Yje 0J e i{f(x e ) tends to the integral II((/0 = J i{/(x)Tl(dx) when the 
sample size (i.e., the cardinality of the index set E) tends to infinity. Let us as- 
sume that the target n has a density n(-) with respect to a reference measure dx. 
Among the Monte Carlo methods, Importance Sampling (see Hesterberg, 1988, 
1995; Ripley, 1987) consists in drawing x e from a proposal probability Q with 
density q(-) and then in weighting this draw with co e = n(x e )/q(x e ) to take into 
account the discordance between the sampling distribution and the target. The 
precision of the importance sampling approximation becomes very unsatisfactory 
when the importance weight oj e has high variance, see Owen and Zhou (2000). 
Such problematic cases arise when the proposal distribution Q is not well adapted 
to the target IT or when the space is of high dimensionality. Practically, the 
fine tuning of the proposal measure for a given target is difficult. However, there 
are adaptive techniques to learn the proposal probability sequentially (see Liu, 
2008; Rubinstein and Kroese, 2004; Pennanen and Koivu, 2004). Among these 
adaptive methods, the PMC scheme of Cappe et al. (2004), generalised in the D- 
kernel paradigm (see Douc et al., 2007a,b; Cappe et al., 2008), aims at building a 
mixture distribution by minimising some optimality criterion (Kullback-Leibler, 
variance,. . . ). Given a parametric family of distributions Q(8) indexed by 6 e 0, 
Cappe et al. (2008) seek the best proposal Q(9*) for a given target II by estimating 
6* sequentially on successive samples. 

In many real problems where computing n(x e ) (hence the importance weight) 
is time consuming, recycling the successive samples generated during the learn- 
ing process is essential. To this end, Cornuet et al. (2012) introduce the Adaptive 
Multiple Importance Sampling (AMIS), combining multiple importance sampling 
methods and adaptive techniques. The AMIS is a sequential scheme in the same 
vein as Cappe et al. (2008). During the learning process, the AMIS tries suc- 
cessive proposal distributions, say Q{9\), . . . , Q(8t)- Each stage of the iterative 
process estimates a better proposal Q(6 t+ i), by minimising a criterion such as, for 
instance, the Kullback-Leibler divergence between Q(9) and II. But the novelty 
of the AMIS is the following recycling procedure of all past simulations. At iter- 
ation t, the AMIS has already produced t samples: x\,...,x l Ni from a first guess 

Q{6\), x\,..., x 2 Ni , from Q(9 2 ), ■ . ., and x\,..., x* N from Q(9 t ) and try to derive a 
new parameter 8 t+ i from all those past simulations. To that purpose, the weight of 
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x] (k < t,i < N k ) is updated with 



.1 * 



(i) 



where 9\,...,6 t are the parameters generated throughout the t past iterations and 
Q, = N\ + N 2 + ■ ■ ■ + N t is the total number of past particles. The importance 
weight (1) is inspired by the techniques of Veach and Guibas (1995). Owen and 
Zhou (2000) popularise those techniques to merge several importance samples 
and even propose a more refined and stabilising alternative, named determinis- 
tic multiple mixture. On various numerical experiments where the target is the 
posterior distribution of some population genetics data sets, Cornuet et al. (2012) 
and Siren et al. (2010) show considerable improvements of the AMIS in Effective 
Sampling Size (denoted further ESS, see Liu, 2008, chapter 2). In such settings 
where calculating n{x k t ) is of high computational cost, the recycling process is the 
appropriate scheme. However, no proof of convergence had yet been provided in 
Cornuet et al. (2012). It is worth noting that the weight (1) introduces long mem- 
ory dependence between the samples, and even a bias which was not controlled 
by theoretical results. The main purpose of this paper is to fill in this gap, and to 
prove the consistency of the algorithm at the cost of a slight modification in the 
adaptive process. Indeed, to get rid of the long memory dependency structure of 
the AMIS, at each iteration t, we suggest learning the new parameter 6 t on the last 
sample x[, . . . , x* Nt weighted with the classical 7r(;c[)|#f jcj, for all i = 1, . . . , N t . 
The only recycling procedure is in the final output that merges all the previously 
generated samples using (1). 

The paper is organised as follows. Our modifications of the original AMIS are 
detailed in Section 2. Our main results are in Section 3, as well as a discussion 
on the assumptions we introduce. The proofs of our results are in Sections 4 and 
5. They rely on limit theorems on triangular arrays (see Chopin, 2004; Douc and 
Moulines, 2008; Cappe et al., 2005, Chapter 9). Multivariate version of those 
limit theorems are given in the Appendix. At last, in Section 6, we compare the 
performance of our new algorithm against the original AMIS and a scheme wit a 
naive recycling strategy. The paper ends with a discussion in Section 7. 
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Algorithm 1 Modified AMIS 

Require: an initial parameter Q\ and increasing sample sizes Ni,. . ., Nr. 

1: for? = 1 T do 
2: for i = 1 — > iVf do 
3: draw x\ from Q(6 t ) 

4: compute co\ = nix'^/qix], 9 t ). 

5: end for 

6: compute 9 t+ i = N; 1 £- =1 ^(^) 

7: end for 

8: set fi T = TVi + • • • + 

9: for ? = 1 T do 
10: for i = 1 — > 7Y ? do 

11: update w< = 7r(^)/Q^ E[=i Wtf 30 

12: end for 

13: end for 

14: return the weighted particles (*}, to}), . . . , (jc;^ , w^), . . . , (x[, to]), (x t Nt ,oj Nt ) 



2 Modifications of the AMIS scheme 

When compared with the recursive algorithm of Cappe et al. (2008), the novelty 
in the AMIS of Cornuet et al. (2012) lies in the recycling process at each itera- 
tion and in the final system which includes all particles generated during the T 
iterations of the adaptive process. Hence, the AMIS estimation of the integral 
II((A) = fi//(x)Tl(dx) is 

" r t=l i=\ 

where £l T = N\-\ h N T is the total number of particles that have been produced 

and 

T 

a/ i = n(x t i )/a- T l J]N k q(x t i ,e k ). 
k=i 

The scheme we propose relies on the same ideas: we base our fit of the pro- 
posal on the same notion of entropy minimisation at each iteration. We assume 
that the minimum of the Kullback-Leibler divergence between the parametric fam- 
ily Q(&) an d the target distribution II lies at 0*. For technical reasons, we further 
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assume that we are able to write 9* as an integral of a known function h(x) with 
respect to the target distribution, namely 9* = f h(x)U(dx). Our process ends with 
the recycling of all past particles by updating all weights with (1). But, contrary 
to Cornuet et al. (2012), the calibration of the new parameter 9 t+1 uses only the 
current sample x\, . . . , x' Nr We hope improvements in the accuracy of 9 t +\ against 
previous estimations by imposing that the sample size N, grows at each iteration. 
Moreover, the influence of the first estimations 9\, . . . , 9 t -\ decreases significantly 
in the mixture of the denominator of the importance weight (1). 

Our modified AMIS is described in Algorithm 1. The learning process be- 
tween lines 1 and 7 draws a sequence of samples from which it calibrates gradu- 
ally the parameter 9. The new value of the proposal's parameter we compute at 
line 6 depends only on the last sample we have drawn. This is the only discrepancy 
from the original algorithm of Cornuet et al. (2012). Here, the recycling process 
producing the final output is silently done by updating the weights of all produced 
samples between lines 9 and 13. Finally, we should note that, if calculating n(x) is 
of high computational cost, the value computed at line 4 can be stored in memory 
to be reused at line 1 1 during the recycling process. 

3 Convergence results 

We state here our main results (on the learning process in Paragraph 3.1, and on the 
final output in Paragraph 3.2). For the sake of clarity, their proofs are postponed to 
Sections 4 and 5. Our set of assumptions is quite complex and so, is discussed in 
Paragraph 3.3, where we provides a much more simple set of assumptions, though 
more restrictive, and a toy example. 

3.1 Convergence of the learning process 

Let us begin with some hypothesis on the parametric family of distributions Q(9) 
and on h(x) which defines the best parameter as 9* = f h(x)U(dx). We assume 
that the space of parameters is a subset of the Euclidean space R d endowed with 
the norm || • ||. From now on, S£ is also a subset of some Euclidean space and 
the reference measure is Lebesgue. We denote q(-,9) : S£ — > BL + the density of 
Q(9) with respect to this reference measure, while n(x) is the density of the target. 
We assume also from now on that IT is absolutely continuous with respect to all 
proposals: V0 e 0, II «: Q(9), which is the minimal hypothesis for importance 



5 



sampling schemes. 

We focus here on the learnt 9 t of our algorithm in order to show convergence 
to the unique optimal parameter 9*. We recall that optimality is defined here in 
terms of Kullback-Leibler divergence. We have already assumed the following on 
9*. 

(KL) The minimum of the Kullback-Leibler divergence between Q(9) and Ti is 
unique, and denoted 

9* = arg min f log I 1 n(x)dx. 

ee& J [q(x,9)) 

Moreover, we have an explicit function h : i?f — » such that 

e-= f« X M^, aM f < + co. 
J J q(x,9*) 

Let us now add some regularity assumptions on the parametric family Q(9), which 
depend on the target n and the function h. 



(Al) On any compact set K c we have 

0eK J [ q{x,9) 

(Al) The integral J n 2 (x)\\h(x)\\ 2 /q(x, 9)dx is finite, and for any n > 0, the map 



r 7r(x)||^(x)|| 

lim sup < — — > n > ax = 0. 



r n 2 {x)\\h{x)\\\ {n{x)\\h{x) u 
9 i — ► I — — — 1^ — - — — — >n}dx 



* 2 (x)\\Kx)\\ 2 . 

q(x, 9) [ q(x, 9) 
is continuous on the the parametric space 0. 

Our first result is a weak law of large numbers in the learning process. 

Theorem 1. Under (Al), the estimate 9 t tends to 9* in probability when t —* oo. 

This first key result allows us to prove the following central limit theorem on 
the learning process, adding assumption (A2). The limit covariance matrix is E, 
given by 



T^ix^hiix^hiix) 

dx-9 9 



q(x, 9*) 
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Theorem 2. Under (Al) and (A2), the scaled discrepancy yN t (9 t -6*) converges 
in law to the multivariate Gaussian distribution ,jV(§, Z) when t — » oo. 

By the way, with this theorem and the Borel-Cantelli lemma, we obtain the 
following strong law of large number, which concludes our results on the learning 
scheme. 

Corollary 3. //'(Al) and (A2) hold and if, moreover, YnT\ 1/N t < +oo, then the 
estimate 9 t converges almost surely to the optimal parameter 9* when t — » oo. 



3.2 Convergence of the final recycling scheme 

The remaining part of our results deals with the final output that merges all the 
samples. More precisely, Theorem 4 below says that the empirical sum of a func- 
tion iff on the merged, re-weighted sample provides a convergent approximation 
of the integral il(i/r). The class of integrands ij/ e L^IT) for which the above 
holds is determined by the following additional assumptions. For any function 
g: — > R, we set \\g\U = sup{|g(*)| : x e 



(A3) The integral f -^M— -r-^ cbc is finite 
J q{x, 9*) 

and the mapping 9 h> 11(G) = \ — — } q{x, 9)dx is continuous over 0. 

¥ J q(x,9*) 



(A4) For any 9 e 0, the integral f { ^(x )/q \x, **** "finite, and 

= J \^(x)il?(x)lq{x,G*)^dx-(Yl(ilf)f. 



we set a, 



(A5) For any rj > 0, the mapping 

6 ^ f -^^q(x,9)l{n(x)\ifr(x)\>T 1 q(x,G*)}dx, 
J q (x, u ) 

"(•¥(•) 



is continuous. 
(A6) The mapping : 9 h-> R^{9) 



is continuous. 



q(;9)q(-,9*) 

We also add the following continuity assumption on the parametric family Q(9). 
(A7) The mapping G, given by G{9) = \\q(-, 9) - q(-, 0*)||oo, is continuous. 
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With all those assumptions at hand, we are now in a position to state conver- 
gence of the final output, which is a strong law of large number. 

Theorem 4. Assume that the parametric family Q(9) satisfies assumptions (Al), 
(A2) and (A7). Assume that c € is a class of functions if/ e h 1 (IT) for which as- 
sumptions (A3) - (A6) hold true. For all the sum over the final weighted 
sample TT^ MIS (if/) given in (2) tends almost surely to J if/(x)n(x)dx when T — » oo. 

3.3 Towards a more simple set of assumptions 

The seven assumptions given below are quite cumbersome to check in concrete 
settings. Besides it is not even obvious that there does not exist any contradiction 
between the seven assumptions (in which case our results would apply only to 
nonexistent cases). The assumption (KL) was already present in the original work 
of Cornuet et al. (2012) and will not be discussed here. 
Let us consider the three following assumptions. 

(Bl) The function h(x) of assumption (KL) is dominated by a polynomial in x, 
i.e., there exists a polynomial P(x) such that \\h(x)\\ < P(x). 

(B2) The density n(x) of the target is a rapidly decreasing function (also called a 
Schwartz function) when \\x\\ — > oo, i.e., for all polynomial P(x), n(x)P(x) — > 
when \\x\\ — » oo. 

(B3) The density q(x, 6) is a function of class ^ 2 in both variables (x, 6) and 
there exists a function P(x, 9), polynomial in x, and whose coefficients are 
continuous in 6, such that 

q(x,8)>l/P(x,6). 

Clearly, assumptions (Bl) - (B3) are much easier to check than (Al) - (A7). 
Assumptions (Al) - (A7) were actually the weakest set of assumptions under 
which we were able to prove our results. But we prove in Lemmas 5 and 6 below 
that the assumptions (Bl) - (B3), together with the fact that the integrand iff is a 
polynomial, imply (Al) - (A7). Certainly, there exists other situations in which 
(Al) - (A7) hold true. Furthermore, Lemma 7 below gives a toy example on 
which (KL), and (Bl) - (B3) hold true. 

Let us briefly discuss (B2) and (B3). It is well known that the basic importance 
sampling scheme is inefficient when the target IT has larger tails than the proposal 
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distribution Q. Indeed, in that case, the weight a>(x) = n(x)/q(x) has high variance 
because q(x) can be much smaller than n(x). The learning process of our algorithm 
begins with a distribution Q(6i), where 8\ is a first guess. Over iterations, this 
process learn gradually a better proposal. Since we do not have much information 
on 6* at the beginning of the process, we might run completely off the road if, by 
accident, we fall on a proposal Q(6) that has smaller tails than IT. Assumptions 
(B3) and (B2) above, controlling the tails of those distributions, ensure that this 
will never happen, and so, are rather moral. Indeed, with (B2), IT has tails that 
decrease much faster than any polynomial while, with (B3), tails of the proposals 
Q(6) cannot decrease faster than polynomially. At last, we shall note that the 
integral J \P(x)\n(x)dx is finite for all polynomials P(x) if (B2) holds and that 
(B2) is fulfilled whenever IT has finite exponential moments. 

Lemma 5. Assumptions (Bl), (B2) and (B3) imply that (Al), (A2) and (A7) hold 
true. 

Proof. Consider a set of random vectors X e , 6 e 0, such that the marginal distri- 
bution of X e is Q(9). Actually, (Al) is nothing but the uniform integrability of the 
set H K = {^ e , 6 e K} of random variables 



when the parameter space is restricted to a compact set K. Using the theorem of 
de la Vallee-Poussin, it is enough to prove that 



since g(x) = x 2 is an increasing function on [0; oo) such that lim^^ g(x) / x = oo. 
But, for any 6, 



is finite because h 2 (x)/q(x,6) is bounded by some polynomial in x and n is a 
Schwartz function. Moreover, using continuity under the integral sign, (4) is a 
continuous function of 9, hence bounded on the compact set K. This proves (3), 
hence (Al). 

Proving (A2) is more technical. Actually, it is enough to prove that, for any 
rj > 0, the function F defined as 



6 = 



\h(X e )\n(X e ) 
q(Xe,6) 




(3) 




(4) 



9 



F(9,6')= f g(x,9)dx, v/i±g(x,e) = 



Ja{6>) 



x\x)\\h(x)\\ 2 
q(x, 8) 



and A{9') = {x : g(x,6') > rj}, is continuous in 6 and 8'. The continuity of 8 h-> 
F(6>, 0') is a classical application of continuity theorem under the integral sign. 
Fix some 8q in the parameter space, and a compact neighbourhood K$ of 8q. With 
(B3), we can bound l/q(x,8) by a polynomial in x. Since the coefficients of this 
polynomial depend continuously on 8, they all can be bounded by constants on 
the compact set Kq. Thus there exists a polynomial P(x) such that, for any 8 e Kq, 
g(x,8) < n 2 (x)\\h(x)\\ 2 P(x) which in integrable because of (Bl) and (B2). This is 
enough to get continuity of 8 i-» F(8, 8'). 

It remains to prove continuity of 8' h-» F(8, at any point 8' of the parameter 
space. We have 



where A(9')AA(9' ) = {x : g(x,ff) > rj > g(x,8' Q )} U {x : gO,^) > r] > g{x,ff)}. 
But, when 8' — » this symmetric difference A(0')AA(S ( ' ) ) tends to the set {x : 
g(x,8 ) = rj}. Indeed, the second difference {x : g(x,8' Q ) > rj > g{x,8')} tends to 
the empty set, since, if 8' — > ff Q , g(x, 8') —* g(x, 8' ) > 77 which is impossible if all 
8' satisfy g(x, 8') < 77. Likewise, when 8' — » O5 the first difference {x : g(x, 8') > 
rj > g(x, d' Q )} tends to {x : O ) = Reporting in (5) leads to 



and this last integral is null because the Lebesgue measure of {x : g(x, 8 ) = t]} is 
null. Which concludes the proof of (A2). 

To prove (A7), we have to prove that G is continuous at any point 8q of the 
parameter space. At least on some neighbourhood of 8q, say K , we have 



where 8 h-» x{8) is continuous. Indeed, x{8) is a solution of d x q(x, 8) - d x q(x, 8*) = 
whose first hand side depends continuously on (x, 6). Hence, using the implicit 
function theorem, its root x{8) is continuous on Using the chain rule in (6), 




(5) 



lim sup \F(8, 8') - F(8, 8' )\ < g(x, 8)l{g(x, = tfdx. 



G(8) = \\q(; 8) - q(; 8*)^ = \q(x(8), 8) - q(x(8), 8*)\ 



(6) 



we obtain continuity of G, thus (A7). 



□ 
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The following lemma gives a class of integrands iff that satisfy our assumptions 
(A3)-(A6). 



Lemma 6. Assume (Bl) - (B3). For any polynomial function ij/(x), assumptions 
(A3) - (A6) hold true. 

Proof. Let us recall that, because of (B2), for any polynomial P{x), J \P(x)\n(x)dx 
is finite. Assumption (B3) and the fact that i/r is a polynomial function implies 
that the integrals of (A3), (A4) and (A5) are all finite. Continuity of the mappings 
given in (A3) and (A5) follows from (B3) and continuity under the integral sign. 
To obtain the continuity under the integral sign at some point 6q, it is enough to 
bound the integrand locally around 9q by an integrable function that not depend 
on 9. Let us take for example the function I*, of (A3). Using (B3), the density 
q(x, 9) might be bounded by some M e in [0; oo) that depends continuously on 9: 

sup q(x, 9) <M e . 

And, by continuity of M g , those constants might be bounded by some M e [0; oo) 
when 9 ranges a compact neighbourhood, say K , of 9 . Thus, the integrand in I*. 
satisfies 

n{x)\fr(x) 7r(x)i/f(x) 
— — — q(x, 9) < M 
q(x,9*) q(x,9*) 

for all 9 in K . This last bound is independent of 9 and is integrable. Which proves 
(A3). We omit the domination of the integrand in (A5) here, but we can obtain it 
following the same type of arguments. Hence (A3), (A4) and (A5) hold true. 

The proof of (A6) follows mainly the same line than the proof of (A7) in 
Lemma 5, and thus is omitted here. □ 

To conclude this discussion on the assumptions, we give below a simple toy 
example where all the above assumptions are true. 

Lemma 7. Assume that the target density IT is the univariate Gaussian distribu- 
tion ,sV(m, 1 ), where m is unknown. Assume also that the parametric family of 
proposals is the set ofStudent-t distribution with 3 degrees of freedom centred on 
9, i.e., 



q(x, 9) = C 



3 



with C = (r(^]V3/r 



Then, h(x) = x and the three simple assumptions (Bl) - (B3) hold true. 
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Proof. Here, n(x) decreases exponentially fast to at infinity, hence (B2). More- 
over, (B3) is trivial because q(x, 9) is already of the form 1 /P(x, 9) where P(x, 9) 
is a polynomial in both x and 9. 

It remains to prove (Bl). Minimising in 9 the Kullback Leibler divergence 
between II and the proposal Q{9) is equivalent to maximising in 9 

J log q(x, 9)n(x)dx = A + Bx A(m - 9), 

where A and B does not depend on 9 and A{u) = Jlog(l + y 2 /3) e~ (y ~ u '?t 2 dy. 
A careful study of the function A(u) shows that it admits a global maximum at 
u = 0. Thus the Kullback Leibler divergence admits a global minimum at 9 = 
m = J xn(x)dx. Hence, (Bl) holds true with h(x) = x. □ 

4 Proofs of the convergence during the learning pro- 
cess 

This Section is devoted to the proofs of Theorems 1, 2 and Corollary 3 which 
are written in Paragraphs 4.1, 4.2 and 4.3 respectively. They rely on the limit 
theorems given in the Appendix. From now on, to distinguish random vectors 
from deterministic vectors, we use the notation Xj instead x\ and we define the 
cr-fields T t = cr{X\ , . . . , X^, . . . , X[~\ ...,X'~ t \) which form a filtration. 

4.1 Proof of Theorem 1 

Sketch of proof. The desired result follows directly from Theorem 14 on the 
triangular array given by 

= n{X\)h{Xj) 
''' N t q{X\A)' 

Indeed, we have 9 t+ \ = £ Ji V t ,i The main difficulty here is in checking assumption 
(Hi) of that theorem. Hence, we begin by proving that 9, is tight, then we check 
that assumption ( Hi) of Theorem 14 is true, and conclude our proof. 

Tightness of 9 t . We have 
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Moreover 



MM 



r, 



= J tt(jc)|| 



h(x) \\dx, 



therefore E(||? (+ i||) < n(||/i||). Then Markov's inequality leads to 

x n (||/z||) 

lim sup P (||0,|| > c) < lim = o. 



Hence, the sequence \6 t } t >\ is tight. 



Checking condition (Hi) of Theorem 14. We have to show that, for any r] > 0, 

the following S , tends to in probability, with 



N, 



s t = J]E[\\v4i{\\v t i>i 1 }\r t \ 



= Jn(4\h(4l\^!^> Nl , l } i , 



q(x, 0,) 



= I(N t Tj,e t ) 



by setting 



Kn,o) 



= J n(x)\\h(x)\\\{ 



{x)\\h{x)\\\\n{x)\Mx)\\>r 1 q(x,e)\&x. 



Fix any s > 0. Using tightness of 6 t proven above, we might introduce a 
compact set K E such that, for all t > 1, F(0 f £ AT £ ) < e. On the event {0, £ AT £ }, we 
have I(N,r],l) t ) < 7(0, 0~ f ) = n(\h\). Furthermore, on the event {6 t e K s }, 



? f 



7(7^, 0,)< sup n(x)\\h(x)\\l 



> N,tj > dx. 



(7) 



Assumption (Al) implies that the right-hand side of (7) goes to since N t rj tends 
to infinity. Thus, 
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limsupE(/(iV f /7,0 t )) < 7r(|/z|)limsupP (0, $ K e ) < n(\h\)e. 

Since e was arbitrary, we have proven that S t — » in L J (P) hence in probability, 
and condition ( Hi) is fulfilled. 

Conclusion. We are now in a position to apply Theorem 14. Condition (i) is 
satisfied because, by construction, given f t , the random variables V t j, for 1 < i < 
N t are conditionally independent and E[||Vy| | T t ] < +oo. Condition (ii) is easy to 
check because, here, 



is not random (thus is a tight sequence). And condition (in) has been proven 
above. Hence the conclusion of Theorem 14 is true and the proof is completed. 
□ 



Theorem 2 follows from a direct application of Theorem 15 on the triangular array 



provided that we can check its four assumptions. 

Checking conditions (i) and (ii) of Theorem 15. By construction, given T t , the 
random variables V t j for 1 < i < N t are conditionally independent by construction. 
Moreover, 




4.2 Proof of Theorem 2 



of 



V t ,i = 



ylN t q{X\A)' 




which is finite by (A2). Hence (i) and (ii) hold true. 
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Checking condition (Hi) of Theorem 15. Fix b e R d , and denote V*., the k th 
component of the vector V t j. On the first hand, we have 



B[<*, V tJ ) 2 \r t ] = £ b 2 k E[ (Vl) 2 \r] + 2 £ M,E[vX]!F r ]. 

fe=l l<k<(<d 

Furthermore, by the continuity assumption in (A2) and Theorem 1, for any pair 
(V*,., V/ f ) of components of V fj! -, we have 



i=i J q\x,9 t J 



J T?{x)h k (x)h e (x) 



Ax in probability. 



q(x, 6*) 
And, on the other hand, 

|](E[<fc.v«>|yv]) a = 2*fl£) a + 2 E 

i=l fc=l \<k<l<d 

Thus, 



2] E[(Z7, y^) 2 |!F f ] - (E[<&, V M >|y>]) — > b'lb in probability, 

i=i 

which proves ( Hi). 

Checking condition (iv) of Theorem 15. Fix a vector e of the canonical basis 
of M, d and a 77 > 0. We have 



Z V ' ml/ w \2-f|/ T7 v| llarl ^ 7T 2 (jc) | | 2 f ^(jc)||/?(x) 



Hence, it is enough to prove that M( y/W t r], 8 t ) tends to in probability, where we 
have set 
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Because of (A2), the function n 2 (x)\\h(x)\\ 2 /q(x, #*)is integrable, hence 

lim M(j], 0*) = 0. (8) 

Using continuity of the map M(rj, •) given in (A2), and Theorem 1 proven above, 

M (i7,? ( ) -> M (r], 6*) in probability. (9) 

Now, fix some s > 0. Because of (8), there exists rjo > such that M (rj , 6*) < s. 
And for all t large enough, yfW t r] > rj thus M( yfW t r], t ) < M(rj , 6 t ). Hence, we 
obtain 

P [M( yfNtTj^t) > 2s] < P [M(77o,?,) > 2s] 

<p[M(77 ,?,)-M(7/ ,r)> e ]. 

This last probability goes to because of (9), whence the convergence in proba- 
bility of M( y/W t r], 6 t ) to for any r] > and ( iv) holds true. □ 

4.3 Proof of Corollary 3 

Because of the convergence in law given in Theorem 2, the random variables 
A ; = N, ||?, - 0*\\ 2 are bounded in L 1 by some finite M , i.e., E\\9 t - 9*\\ 2 < M /N t . 
Then, with Chebyshev's inequality, this yields 

F(\\e t -e*\\>s)<M (sN t y l . 

The Borel-Cantelli lemma leads to the almost sure convergence of 6, because 
YjtMoisNtY 1 is finite by assumption. 

5 Proof of the convergence of the final recycling scheme 

The last step is to prove Theorem 4, i.e., that the AMIS estimator 

n£ MIS W = ^ f f n(X '\ 
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which is the result of the unique recycling step of our scheme, is consistent for the 
integral il(</0 for a large class of functions iff. 

The proof is organised as follows. We introduce a sequence of random vari- 
ables, Il^i/O and prove that this sequence tends to H(i//) almost surely. We then 
show that the difference between our estimator and this auxiliary variable, namely 
n^ MIS (i/0 - U* T (tf/), tends to almost surely. Hence the consistency stated in The- 
orem 4. 

5.1 Convergence of some auxiliary variable 

Let us introduce the sequence 



(which cannot be computed in practice because 8* is unknown). And note that this 
auxiliary variable is a (weighted) average of the random variables 7rJ(^) given by 



Proposition 8. Assume (Al) - (A5). When t — » go, ^fW t {n* t (tff) - /!(0 f _i)} tends in 
law to the Gaussian distribution ^V(0, <x^), where 1^(6) is the mapping defined in 
Assumption (A3). 

Proof. The approach here is very similar to the proof of Theorem 2. Indeed, the 
desired result follows from a direct application of Theorem 15 on the triangular 
array given by 



Likewise, the only challenge is in checking the four conditions in order to apply 
this theorem. 

By construction of the Xj, the V t j are independent conditionally on T t - Thus, 
condition ( i) is satisfied. Moreover, 





Let us begin with a result on this last sequence of random variables. 



V t ,i = 



7T(g)jKg) 
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and this last integral is finite because of (A4). Thus, condition (ii) is fulfilled. 
Using Theorem 1 as well as continuity arguments on the mappings of (A3) and 
(A5), 

|{E[|y,| 2 |^-(E[y,|^) 2 }= / Jg^^ 

-{J -^¥) qM)dx ) 

— > cr^ in probability. 

Which proves condition (Hi). 

It remains to prove that condition (iv) is also fulfilled. For any fixed rj > 0, we 
have 



Thus, in order to prove condition (iv), it is enough to show that 

M*(y/N t r],'e t ) — »0 in probability, (10) 

where we have set 

Because of (A3), the function n(x)\tfr(x)\/q(x, 6*) is integrable which implies that 

limM*O7,0*) = O (11) 

Tj— >OC> 

Besides, using continuity of the mapping given in (A3), and the convergence of 6, 
proven in Theorem 1 , we obtain 

AT(77,?,) — > M*(7],6*) in probability. (12) 

Now fix s > 0. Using (11), there exists i] > such that M*(rj Q , 6) < s. And 
for all t large enough, yfW t r] > r/ , thus 
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We have then 

P [M*( y/NtriA) >2s]<¥ [M*(i]oA) >2s]<¥ [M*( m ,9 t ) - M*( m , 9*) > s] 

and this last probability goes to because of (12). Since s is arbitrary, this leads 
to 

M*( ^jN t flA) — > in probability. 
Whence (10), condition (iv) and the desired result by applying Theorem 15. □ 

We shall now recall that the auxiliary variable H* T (iff) is a weighted average of 
the n^(ifr) for t = 1, . . . , T which are controlled by Proposition 8 proven above. 
Thus, the following generalised Cesaro Lemma will be useful. We omit the proof 
which is classical. However, we should take care of the fact that there exists 
no Cesaro Lemma for the convergence in probabilty, see for instance Billingsley 
(1995), p. 543 and following. 

Lemma 9. If{U t ) is a sequence of random variables that converges to £/«, almost 
surely and if{b t } is a sequence of positive real numbers such that B t = b\ + . . . + b t 
tends to infinity, then B~ l Y!k=\ bfU t converges to £/«, almost surely. 

This Cesaro Lemma and Proposition 8 above leads to the following. 

Proposition 10. If '(Al) - (A5) hold true and ifj^ 1/N t < +oo, then 

n^((/0 — > n(i/r) almost surely. 
Proof. Proposition 8 implies that 



K := supMEfeOA) - i;(e t )} 2 ] < +oo. 

f>i 



£ p flaw - >s)< e~ 2 Kj] ^ < 



And for any s > 0, Tchebychev inequality yields 

+00 

Applying the Borel-Cantelli lemma, it follows that n* t (\]/) - I^(6 t ) —> almost 
surely. Using (A3) and Corollary 3, Il(6 t ) tends to n(i/^) almost surely, thus n* t (ifr) 
too. Applying Lemma 9 with b t = N t yields the convergence of TT T (iff). □ 
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5.2 Controlling the discrepency between the AMIS estimator 
and the auxiliary variable 

The convergence of Hj MIS (iff)-U* T (if/) towards almost surely is in Proposition 12 
below, whose proof relies on some preliminary result given in Lemma 1 1 . To this 
end, we define the function D(-, 6\ j) : 3£ i-» K.+ by 



D(x,9 UT ) = cV j]N k q(x,e k ) 



k=l 

which appears in the denominator of the updated weight (1). Because of the con- 
sistency of the learning scheme proven in Section 4, we are able to show in the 
following lemma that this denominator resembles the denominator of the classical 
importance weight, when the proposal distribution is Q(9*). 

Lemma 11. If assumptions (A.l) (A2) and (A7) are fulfilled, and if moreover, 
2 f =7 l/N t < +oo, then 



lim 



D(-,9i;t) - q (-,&*) =0 almost surely. 



Proof. We have 



D(;9 1:T )-q(-,9*) 



<Q^ Y,N t \\q(;9 t )-q(;e*) 

t=\ 

T 

1 

— 1 



(13) 



Because of (A7), G is continuous. Thus, using Corollary 3, G(9 t ) —* almost 
surely. Applying Lemma 9 with U t = G(9 t ) and b, = N t on the right hand side of 
(13) leads to the desired convergence of the distance \\D(-, 9\,t) - <?(•, ^*)IU. □ 

We can now state and prove the result controlling the difference between the 
AMIS estimator and the auxiliary variable T\* T (if/). 

Proposition 12. If conditions (Al), (A2), (A6) and (A7) hold true, and if, more- 
over, 2,=°i 1/N t < +oo, then 



lim nf" 5 (</0 - Tl* T (if/) = almost surely. 
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Proof. It follows from Corollary 3 and from the continuity of the function 
given in (A6) that R^(6 t ) — » R^{9*) almost surely. With Lemma 9, this yields 
Q,j l Yik=i NkR^{6k) —* R<ji{&*) almost surely. Hence, with Lemma 11, 



lim 

T— *+oo 



V k=\ 



D(;e UT )-q(;9*) = almost surely. (14) 



On the other hand, we have 



nf^W-rLM-A) 



T N, 



< 



< 



Lj f=l i=\ 

! T N, 

Lee 



D{x\,e x , T ) q(%,0*) 

7r(Xf)lfr(Xf) 



D{X\,e hT )q(4,6*) 



\D(xj,e hT )-q(xl,e*) 



The difference D(X ( ',0i :r ) - g(X},0*) might be bounded by | D(-,6i :T ) - q(-,0*) 
Moverover, convexity of the function ih 1/x on R+ proves that the ratio in each 
term of the above sum might be bounded by 



D(x\;e\. T )q{x\,e*) < 



Therefore 



|nr s W-n* r w|<^-X£ £ 



w < ( L. N k RM) 



D(;e 1:T )-q(-,e*) 



t =i i=i yk=i 

Then the right-hand side of this display tends to almost surely because of Lemma 9 
and (14). Which concludes the proof of Proposition 12. □ 



5.3 Conclusion of the proof of Theorem 4 

Proposition 10 gives the convergence of the auxiliary variable H* T (if/) towards the 
integral II(i/0 almost surely, while the discrepency between the AMIS estimator 
and this auxiliary variable becomes negligeable almost surely (Propositon 12). 
Whence the almost sure convergence of the estimator IIy MIS (i/r) to the integral 
n(t/0, and then the proof of Theorem 4 is completed. □ 
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6 Numerical experiments 



We provide here two numerical examples on which we have compared 

[a] an iterative and adaptive algorithm learning 6 with a naive recycling strategy 
at the end; 

[b] the original AMIS of Cornuet et al. (2012) and 

[c] our modified AMIS with its recycling step only at the end. 

We refer the reader to Cornuet et al. (2012) for the orginal algorithm [b]. The 
algorithm [c] is the main topic of this paper, and is described with great details in 
Section 2. At last, algorithm [a] follows the same code lines, but stops at line 7 
and returns a merging of all past samples without updating the weights computed 
at line 4. 

The first target, which is described in Paragraph 6.1, is the banana-shaped tar- 
get introduced by Haario et al. (1999), who is known to be a difficult target. In the 
second example, which is in Paragraph 6.2, the target is the posterior distribution 
when conducting a Bayesian analysis on a population genetic data set. 

In both cases, it turns out that our algorithm was the most powerful for a 
given computational cost, that is to say, for a given number of simulations from 
proposals. 

6.1 Performance on a benchmark target 

The banana-shaped target of Haario et al. (1999) is based on a multivariate Gaus- 
sian law jV(Q, S), with covariance matrix 2 = diag(<x 2 , 1, . . . , 1), that is non 
linearly transformed on the second component. The non linear transformation 
depends on a parameter b > and defines the target 

n(x) = f(x u x 2 + b(x\ - a 2 ), x d ). 

In this simulation study, we have set b = 0.03 and cr = 10. The corresponding 
density is given in Figure 1 when d = 2. 

The previously described schemes [a], [b] and [c] rely all on a family of pro- 
posals. Because of the shape of the target, we are using mixtures of 6 multivariate 
Student-? distributions as proposals Q(6), whose densities writes 



6 




(15) 
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Figure 1: Numerical results of the three algorithms [a], [b] and [c] on the banana-shaped 
target. Top-left: the density of the strongly distorded banana-shaped target through the 
contour lines of confidence areas with probability 0.5, 0.9 and 0.99. Top-right: values 
of the ESS over 100 replicates of the three algorithms. Bottom-left and Bottom-right: 
boxplots representing the variability (over 100 replicates) when computing the variances 
of the two first marginal distributions of the target with the three algorithms. The ttue 
marginal variances of the target correspond to the horizontal, dotted lines. 
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where t v (x,fi, S) is the f-distribution with 3 degrees of freedom, centered at /i and 
with dispersion matrix E. This family of proposals is then parametrised by 

e = (p u . ..,p k ,nx,.. . ,Mk,^u ■ ■ • 

At each iteration t of all three algorithms, the new parameter # f+1 is obtained by 
fitting the mixture (15) on the weighted sample approximating the target. In algo- 
rithms [a] and [c], this weighted sample is the current sample, while [b] performs 
a recycling on all past samples. And the fitting procedure is the EM scheme of 
Peel and McLachlan (2000). All three algorithms were started with a same first 
guess 6\ computed as in Cornuet et al. (2012). 

Numerical results on one hundred replicates of all three algorithms are re- 
ported in Figure 1 when the number of iterations is T = 36 and the number of 
particles at iteration t is N t = 5,000 + 1,000?, thus a total number of 846,000 
draws from proposals. It turns out that all three algorithms are similar in term of 
error when estimating the variances of the first two marginal distributions of the 
target. But, clearly, the ESS of the merged samples produced by algorithms [b] 
and [c] are much higher than the ESS of algorithm [a] that suffer from the (poten- 
tially) high variances of the un-corrected weights iiix^lqix 1 ^ 6 t ). And there seems 
to be a small benefit for algorithm [c] in terms of ESS when compared with [b]. 

Finally, we have tried other ways to allocate the overall computational cost, 
and for instance, we have tried sample sizes N, growing mainly as t 2 (so that 
our theoretical assumption Yit ^l^t < oo is true). The corresponding numerical 
results are not presented in detail here, but they were very similar to the results of 
Figure 1. 

6.2 A population genetics' simulation study 

Our second example comes from a population genetic problem. More precisely, 
we want to conduct a Bayesian analysis of a genetic data set *2) to infer mutation 
and migration rates in a parametric model. Assume that the species of interest is 
composed of two large populations at equilibrium, one on an island and the other 
one on the mainland. The parametric model we have used is coalescent based, and 
is detailed, for instance in Donnelly and Tavare (1995) and Rousset and Leblois 
(2012). The genetic data come from two samples of individuals corresponding 
respectively to the two populations, genotyped at five independent micro satellite 
loci. 

Hence the target n{x) is the posterior distribution on a two dimensional pa- 
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Figure 2: Comparison of the three algorithms [a], [b] and [c] over 100 replicates based on 
distances between the discrete probability function given by the output and the target in the 
population genetic example. Left: Cramer-von-Mises distances. Middle: L 2 -distances. 
Right: L°°-distances. 

rameter x = (x^g, x mut ) when the prior is a non informative uniform distribution 
on some compact set Y. This example is actually typical of situations where the 
density of the target n(x) is of high computational cost. With coalescent based 
models, the likelihood, thus the posterior density n(x), is an integral over a latent 
process that, by chance, is computed via importance sampling too, see De Iorio 
et al. (2005) and Rousset and Leblois (2012). 

The family of proposals in the three sequential algorithms is composed of bi- 
variate Gaussian distributions, conditioned (or truncated) on the support Y of the 
prior distribution. The parameter of the proposal is a four dimensional vector 
9 = (/Vig,AW, cr^o'mut)' whose first two components give the position of the 
mode and last two components give the marginal variances of the diagonal covari- 
ance matrix. 

Performances of the sampling algorithms were compared as follows. Instead 
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of looking at the estimates of II(t/0 for various integrands iff, we decided here 
to evaluate the outputs with distances between the discrete measure induced by 
the weighted final samples and the target. In Figure 2, we have represented the 
Cramer-von Mises, the L 2 - and the L ro - distances between the empirical distri- 
bution function F T (x) of the final weighted samples and the distribution function 
F(x) of the target, i.e., the posterior distribution. Furthermore, since the density of 
the target cannot be written with a close formula in the concrete example, we shall 
describe how F(x) was computed. Following an idea of R. Leblois and F. Rousset, 
estimates of n(x) were computed for values of x ranging a regular 500 x 500 grid 
of the support of the prior distribution. The estimation error was then decreased 
using a kriging model on n(x), assuming regularity conditions of that posterior 
density. Of course, this sharp approximation comes at a much higher computa- 
tional cost than any run of the Monte Carlo algorithms we compare here. 

The results presented in Figure 2 exhibit a clear advantage to our modified 
AMIS, namely [c], in front of [a], the sequential scheme with a naive recycling 
and [b], the original AMIS, whatever the distance. Thus, modifying the original 
AMIS was not only a way to obtain the theoretical results of Section 3, but also a 
real improvement of the original algorithm. One of the reason that might explain 
this phenomenon is that the recycling scheme of the original AMIS introduces a 
bias on 9 during the learning process which tends to accumulate, and thus is large 
enough to degrade the output quality when compared to our modified AMIS. 

7 Discussion 

For a certain class of functions, we derived a convergence result of our modi- 
fied AMIS. We have proved a strong law of large numbers for a large class of 
integrands characterized by regularity conditions and for a general family of pro- 
posals. In future research, we could try to improve on our results by relaxing these 
hypotheses, using other probabilistic tools to obtain a central limit theorem on the 
final output, or at least a rate of convergence. Our technique of proof fits the case 
where N t tends to oo. In other words, we allocate more computation time at the 
end of the learning process, when a good proposal has been obtained. Which is 
not completly crazy. But £ 1 /N t < oo is a technical assumption we could try to 
remove in the future. 

Besides, another route might be taken to prove theoretical results on the mod- 
ified AMIS, based on Markovian arguments. Indeed, when the sample size does 
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not vary between iterations during the learning process, i.e., when N\ = N% = 
. . . = Nt = N, the sequence of pairs (x\. N ,0 t ) form a Markov chain. And the 
final sum in (2) might be traited using results on averages over a path of a Markov 
chain. But we have left this route for future works. 

The numerical experiments presented above exhibit an advantage to our mod- 
ification of the AMIS at least in two situations. To decide which one is the best 
among the original and the modified AMIS, a more careful comparison on differ- 
ent targets is required. For instance, Cornuet et al. (2012) and Siren et al. (2010) 
gave more accomplished numerical studies with the original AMIS scheme in 
population genetics. We can conduct the comparison study on both examples. 

Finally two important, methodological issues have not been tackle in this the- 
oritical work. The first one deals with the initialization of the AMIS. The original 
paper of Cornuet et al. (2012) proposed an answer based on a logistic sample 
when nothing is known on the target. We stress here that the starting distribution 
is of great practical consequence: for instance, if the first sample misses a mode 
of the target distribution, we have almost no chance to see it during the whole 
process. That was summed up by Cornuet et al. (2012) as the "what-you-get-is- 
what-you-see" nature of the AMIS. Likewise, a recurring numerical question on 
the AMIS concerns the allocation of the overall computational cost (given by the 
final system size Q T ). To optimize allocation, one could propose and study alloca- 
tion strategies on different iterations via the sequence Ni, . . . , N T . We also believe 
that the winner in the competition between the original and the modified AMIS 
depends also on this allocation strategy 

Appendix - Limit results on triangular arrays 

This appendix summarizes two results on the asymptotics of triangular arrays that 
are widely used in our proofs. These results are multidimensional generalisations 
of two theorems given in Chapter 9 of Cappe et al. (2005). In the following, all 
random variables are assumed to be defined on a joint probability space (S, T , P) 
and {N t } t>l denotes an increasing sequence of integers. In those theorem as well 
as throughout this paper, we used the notion of tightness of random variables that 
we recall here. (See Billingsley (1995), p. 336 for more details) 

Definition 13. A sequence of random vectors {£/„} is tight if 
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lim supP(||£/J >7]) = 0. 

7^°° «>1 V ' 



We have the following law of large numbers and central limit theorems. 



Theorem 14 (Law of large numbers). Let yt,\\ x< . <N be a triangular array of 

random vectors on W 1 , and let {T t ) t>l be a sequence of cr -fields. Assume that the 
following conditions hold true. 

( i) The V ti for i = l,...,N t are conditionally independent given Tt, and for any 



v4 



< +00. 



t and i = 1 , . . . , N t , E 

(ii) The sequence j E[||v f , f || T t fj is tight. 
( Hi) For any positive 77, 

N, 



2>[N|i{lNH 



T, 



Then, 



(=1 



N, 



in probability. 



J] |v f>j - E[v M |y,]J — * in probability. 
(=1 

Theorem 15 (Central limit theorem). Let {Vt,i} l<<N be a triangular array of 

random vectors on R d , and let ( < F t ) t>l be a sequence of cr- fields. Assume that the 
following conditions hold true. 

( i) The V ti for i = 1, . . . ,N t are conditionally independent given T t . 

(ii) For any t > 1 and i = l,...,N t , E[||V f ,,|| 2 |!F f ] < +00. 

( Hi) For any vector b e R d , there exists a constant cr 2 h > such that 

J] |e[<£, V tJ ) 2 \r t ] - (E[(b, V tJ )\r t ]) } oj, in probability. 
(=1 
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( iv) For all vector e of the canonical basis ofR d and any positive number n 



converge in distribution to ,yV(0, Z) when t — > +00. 
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